###############################################
## Code for Analysing WTO and Free Trade
## Generate Figures 7, A12, A15, A16, A17 and Table A3
###############################################
rm(list = ls())


sink("log/log_WTOReplication.txt") # start log
begin.time<-Sys.time()


## compile functions 
library(bpNet)
#library(mcmcr)
library(abind)
library(ggplot2)
library(ggstance)
library(stargazer)
library(coda)
library(igraph)


aa <- -0.6  ##??
bb <- 0.1    ## ??
set.seed(123456789)

## ------------------------- ##

## setwd('~/Desktop/PSRMreplicationLP/data/EmpData')
source('code/net_plot.R', chdir = TRUE)    ## load code to make figures
## Load data
traderep <- read.csv("data/EmpData/tradewlag.csv")
## set the missing data as NA rather than -999
traderep$l1yrsoffc[traderep$l1yrsoffc==-999] <- 0
rhoZ1 <- read.csv("data/EmpData/WTOdensity.csv")[,-1]
load("data/EmpData/WmWTOlist.RData") # WTO network

attach(traderep)
syear <- 1971
eyear <- 2008
Tt <- eyear-syear+1
## original scale of tariff
Y0 <- newtar[year>(syear-1)&year<(eyear+1)]
## DV in the paper 
Y <- log(newtar[year>(syear-1)&year<(eyear+1)]+0.01) 
n <- length(Y)
X1 <- cbind(l1polity2, l1log.pop, l1log.gdppc,l1bpc,l1ecris,l1yrsoffc,
              l1gatt.wto, l1usheg, l1avetar, Eight, Wmember, wtrade, lagnewtar)[year>(syear-1)&year<(eyear+1),]
X0 <- cbind(l1polity2, l1log.pop, l1log.gdppc, l1bpc, l1ecris, l1yrsoffc, l1usheg, l1avetar, Eight, wtrade)[year>(syear-1)&year<(eyear+1),]
Uid <- cowcode[year>(syear-1)&year<(eyear+1)]
Tid <- traderep$year[year>(syear-1)&year<(eyear+1)]
## put together the data
dataWTO <- data.frame(Uid, Tid, Y, X1)

tdataWTO <- data.frame(Y0, X0)
names(tdataWTO) <- c("TARIFF", "REGIME", "LN POP", "GDP PC",
                     "BP CRISIS", "EC CRISIS", "OFFICE", "US HEGEMONY", 
                     "AV TARIFF", "8 OPEN", "GLOBAL WTO TRADE")


## produce Table A3
stargazer(tdataWTO)

## produce the final row of Table A3 
tdensity <- data.frame(rhoZ1)
names(tdensity) <- "DENSITY"
stargazer(tdensity)

## -----Clean the network data and generate W matrix ------------ ## 
W <- WmWTOlist[-1]
W2 <- array(as.numeric(unlist(W)), dim=c(85, 85, 38))[1:85,1:85,1:Tt]

## Get dimension of TSCS 
N <- dim(W2)[1]
TT <- dim(W2)[3]

## correct the W matrix (## The matrix has nonzero diagnal elements )
## and standardize the matrix
 W3 <- W2
 for (i in 1:TT){
    w <-  W2[ , , i]
    diag(w) <- 0
    ww <- (w/rowSums(w))
    ww[ww=="NaN"] <- 0
    W3[ , , i] <- ww
    }

## combine the network measure and group-level intercept
rhoZ <- cbind(1, rhoZ1)[1:TT,]
#int <- c(1, 1, rep(0, (TT-2)))[1:TT,] ## two 1's because we include y lag and the first period will be taken out
#rhoZ2 <- cbind(int, rhoZ1)[1:TT,]

zero <- cbind(rep(0, TT))

## number of total iterations and number of burnin
mcmc <- 105000  ## this analysis needs much more iterations to converge
burnin <- 5000
factorn <- 20   # set the number of factors

##-------- Simulation begins ----------##

## rho_t-MLST-MF: the proposed model
outMF <- bpNet(data = dataWTO,   
	              W = W3,         
	              index = c("Uid", "Tid"),
	              Yname = "Y",
	              Xname = c("l1polity2", "l1log.pop", "l1log.gdppc","l1bpc","l1ecris","l1yrsoffc",  "l1gatt.wto"),
	              Zname = c("l1usheg", "l1avetar", "Eight",  "wtrade"),  
	              Aname = NULL,  
                      Contextual = c("l1bpc","l1ecris"), 
	              Contextual.effect = "both",  
	              lagY = TRUE,               
	              force = "both",             
	              r = factorn,          
	              flasso = 1,             
	              rhoZ = as.matrix(rhoZ),      
	              rho_spec = "multilevel",                
                      niter = mcmc,
                      burn = burnin)

##### Model rho_t-MLST-FE: Run the FE model
## Try whether we need to reload the data as in the folder of WTOreplication2022

outFE <- bpNet(data = dataWTO,   
	              W = W3,         
	              index = c("Uid", "Tid"), 
	              Yname = "Y",
	              Xname = c("l1polity2","l1log.pop", "l1log.gdppc","l1bpc","l1ecris","l1yrsoffc", "l1gatt.wto"),
	              Zname = c("l1usheg", "l1avetar", "Eight",  "wtrade"),  
	              Aname = NULL,  
                      Contextual = c("l1bpc","l1ecris"), 
	              Contextual.effect = "both", 
	              lagY = TRUE,              
	              force = "both",           
	              r = 0,          
	              flasso = 0,             
	              rhoZ = as.matrix(rhoZ),       
	              rho_spec = "multilevel",         
                      niter = mcmc,
                      burn = burnin)


## Model rho-MLST-MF: constant network interdependence

outCon <- bpNet(data = dataWTO,  
	              W = W3,         
	              index = c("Uid", "Tid"), 
	              Yname = "Y",
	              Xname = c("l1polity2","l1log.pop", "l1log.gdppc","l1bpc","l1ecris","l1yrsoffc", "l1gatt.wto"),
	              Zname = c("l1usheg", "l1avetar", "Eight",  "wtrade"),   
	              Aname = NULL,  
                      Contextual = c("l1bpc","l1ecris"),  
	              Contextual.effect = "both", 
	              lagY = TRUE,               
	              force = "both",            
	              r = factorn,         
	              flasso = 1,           
	              rhoZ = as.matrix(rhoZ),      
	              rho_spec = "constant",       
                      niter = mcmc,
                      burn = burnin)
##------------------Simulations end here -----##

####################
## Draw graphs
####################

## Figure A15: factor selection
p1 <- omegaPlot(outMF, oxlim = c(-1, 1), plotlength = 5, labelname = "Factor")
ggsave('plots/WTO/factorslog.pdf', p1, width = 18, height = 12)

## Figure A16: factor selection
p1 <- omegaPlot(outCon, oxlim = c(-1, 1), plotlength = 5, labelname = "Factor")
ggsave('plots/WTO/factorslogCon.pdf', p1, width = 18, height = 12)

### Figure 7
### Figure 7 (a): rho
rho <- rep(0, (TT-1))
modelRho <- outMF$Rho
## ----------------------- 1 varying rho_t
datap <- cbind.data.frame((syear+1):eyear, apply(modelRho, 1, mean),
		          apply(modelRho, 1, quantile, 0.025),
		          apply(modelRho, 1, quantile, 0.975),
		          c(rho))
names(datap) <- c("time", "mean", "cil", "ciu", "true")
datap$type <- rep("m", (TT-1))
p <- ggplot(datap) + xlab("Time") + ylab("Rho")
p <- p + geom_line(aes(time, mean, colour = type, size = type, linetype = type))
p <- p + geom_line(aes(time, true), size = 0.5, colour = "red", linetype = "dashed")
p <- p + geom_ribbon(aes(x = time, ymin = cil, 
	              ymax = ciu), fill = "#000000", alpha = 0.2)
set.limits <- c("t", "m", "c")
set.labels <- c("0", "Posterior Mean", "95% CI")
set.colors <- c("red", "grey50", "#00000080")
set.size <- c(1, 1, 3)
set.type <- c("dashed", "dashed", "solid")
p <- p + scale_colour_manual(limits = set.limits,
	                             labels = set.labels,
	                             values =set.colors) +
	         scale_size_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.size) +
	         scale_linetype_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.type) +
	         guides(colour = guide_legend(title=NULL, nrow=1),
	                size = guide_legend(title=NULL, nrow=1),
	                linetype = guide_legend(title=NULL, nrow=1))
p <- p + theme(legend.text = element_text(margin = margin(r = 10, unit = "pt"), size = 12),
	               legend.position = "bottom",
	               axis.title = element_text(size=12),
	               axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0)),
	               axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
	               axis.text = element_text(color="black", size=8),
	               axis.text.x = element_text(size = 8),
	               axis.text.y = element_text(size = 8),
	               plot.title = element_text(size = 15,
	                                         hjust = 0.5,
	                                         face="bold",
	                                         margin = margin(10, 0, 10, 0)))
p <- p + scale_y_continuous(limits = c(aa, bb))+ xlab(" ") + ylab(" ")
a <- 25
p <- p +theme(axis.text.x = element_text( hjust = 1, size=a))
p <- p +theme(axis.text.y = element_text( hjust = 1, size=a))
p <- p+theme(legend.title = element_text(size=25), legend.text = element_text(size=20), legend.key.size = unit(2, 'cm'))
ggsave('plots/WTO/RhoMF.pdf', p, width = 9, height = 8)
## Save another one with different proportions 


## rho_t-MLST-MF
## posterior mean in 1972
datap$mean[1]
## posterior mean in 2008
datap$mean[ (TT-1)]
## the standard deviation of posterior mean
sd(datap$mean)


cat("Model rho_t-MLST-MF: 1972 and 2008 \n")
print(datap$mean[c(1, (TT-1))])
cat("\n")




### Figure 7 (b): rho
rho <- rep(0, (TT-1))
modelRho <- outFE$Rho
## ----------------------- 1 varying rho_t
datap <- cbind.data.frame((syear+1):eyear, apply(modelRho, 1, mean),
		                          apply(modelRho, 1, quantile, 0.025),
		                          apply(modelRho, 1, quantile, 0.975),
		                          c(rho))
names(datap) <- c("time", "mean", "cil", "ciu", "true")
datap$type <- rep("m", (TT-1))
p <- ggplot(datap) + xlab("Time") + ylab("Rho")
p <- p + geom_line(aes(time, mean, colour = type, size = type, linetype = type))
p <- p + geom_line(aes(time, true), size = 0.5, colour = "red", linetype = "dashed")
p <- p + geom_ribbon(aes(x = time, ymin = cil, 
	             ymax = ciu), fill = "#000000", alpha = 0.2)
set.limits <- c("t", "m", "c")
set.labels <- c("0", "Posterior Mean", "95% CI")
set.colors <- c("red", "grey50", "#00000080")
set.size <- c(1, 1, 3)
set.type <- c("dashed", "dashed", "solid")
p <- p + scale_colour_manual(limits = set.limits,
	                             labels = set.labels,
	                             values =set.colors) +
	         scale_size_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.size) +
	         scale_linetype_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.type) +
	         guides(colour = guide_legend(title=NULL, nrow=1),
	                size = guide_legend(title=NULL, nrow=1),
	                linetype = guide_legend(title=NULL, nrow=1))
p <- p + theme(legend.text = element_text(margin = margin(r = 10, unit = "pt"), size = 12),
	               legend.position = "bottom",
	               axis.title = element_text(size=12),
	               axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0)),
	               axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
	               axis.text = element_text(color="black", size=8),
	               axis.text.x = element_text(size = 8),
	               axis.text.y = element_text(size = 8),
	               plot.title = element_text(size = 15,
	                                         hjust = 0.5,
	                                         face="bold",
	                                         margin = margin(10, 0, 10, 0)))
p <- p + scale_y_continuous(limits = c(aa, bb)) + xlab(" ") + ylab(" ")
a <- 25
p <- p +theme(axis.text.x = element_text( hjust = 1, size=a))
p <- p +theme(axis.text.y = element_text( hjust = 1, size=a))
p <- p+theme(legend.title = element_text(size=25), legend.text = element_text(size=20), legend.key.size = unit(2, 'cm'))
ggsave('plots/WTO/RhoFE.pdf', p, width = 9, height = 8)

## rho_t-MLST-FE
## posterior mean in 1972
datap$mean[1]
## posterior mean in 2008
datap$mean[ (TT-1)]
## the standard deviation of posterior mean
sd(datap$mean)


cat("Model rho_t-MLST-FE: 1972 and 2008 \n")
print(datap$mean[c(1, (TT-1))])
cat("\n")

## ### Figure 7 (c): rho
rho <- as.mcmc(t(outCon$rho0))
In <- rep("rho", dim(rho)[1])
rhop <- data.frame(rho, In)
p <- ggplot(data= rhop , aes(x=In, y=rho))+
    geom_boxplot() + xlab(" ") + ylab(" ")
p <- p + scale_y_continuous(limits = c(aa, bb))  
a <- 25
p <- p +theme(axis.text.x = element_text( hjust = 1, size=a))
p <- p +theme(axis.text.y = element_text( hjust = 1, size=a))
#rhoA <- as.mcmc(t(out$rhoA))
ggsave('plots/WTO/RhoCon.pdf', p, width = 3, height = 8)
mean(rho)

cat("Model rho-MLST-MF: constant \n")
print(mean(rho))
cat("\n")



## Figure A17: coefficients in 3 models

paraU1 <- as.mcmc(t(outMF$beta))
paraG1 <- as.mcmc(t(outMF$rhoA))
paray1 <- as.mcmc(t(outMF$rhoy))
paraU2 <- as.mcmc(t(outFE$beta))
paraG2 <- as.mcmc(t(outFE$rhoA))
paray2 <- as.mcmc(t(outFE$rhoy))
paraU3 <- as.mcmc(t(outCon$beta))
paray3 <- as.mcmc(t(outCon$rhoy))
aa <- rbind(summary(paray1)[[2]], summary(paraU1)[[2]], summary(paraG1)[[2]], summary(paray2)[[2]], summary(paraU2)[[2]], summary(paraG2)[[2]], summary(paray3)[[2]], summary(paraU3)[[2]])

model <- c(rep("WTO", 13), rep("FE", 13), rep("Constant", 11))
nameC <-  c("Y_t-1", "Intercept", "REGIME", "LN POP", "GDP PC","BP CRISIS","EC CRISIS","OFFICE",  "GATT",  "W*BP CRISIS","W*EC CRISIS", "rho0", "DENSITY")
nameCCC <-  c("Y_t-1", "Intercept", "REGIME", "LN POP", "GDP PC","BP CRISIS","EC CRISIS","OFFICE",  "GATT", "W*BP CRISIS","W*EC CRISIS")
parameter <- c(nameC, nameC, nameCCC)
ddd <- data.frame(parameter, aa, model)
names(ddd) <- c("parameter", "low1", "low2", "Mean", "Upper2", "Upper1", "model")

## put those posteriors at similar scale in one graph
nn1 <- c(2, 10, 11)

## Figure A17: middle panel
upa <- ddd[ c(nn1,(nn1+13),(nn1+13+13)),]

CIplot1 <- ggplot(upa, aes(colour = factor(model)))
CIplot1 <- CIplot1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2)
CIplot1 <- CIplot1 + geom_linerangeh(aes(y = factor(parameter),
                                         xmin = low2, xmax = Upper2), lwd =1, position = position_dodgev(height = 1/2))

CIplot1 <- CIplot1 + geom_pointrangeh(aes(y = factor(parameter)
                                          , x = Mean, xmin = low1, xmax = Upper1), lwd =.45,  position = position_dodgev(height = 1/2),
                                       fill = "WHITE")
CIplot1 <- CIplot1 + theme_bw()
CIplot1 <- CIplot1 + theme(legend.position="bottom")
CIplot1 <- CIplot1 + labs(colour = "Model")
CIplot1 <- CIplot1+ scale_colour_discrete(labels = c("rho-MF", "rho_t-FE",   "rho_t-MF")) + xlab(expression(beta)) + ylab("")+ggtitle("")
a <- 16
CIplot1 <- CIplot1 +theme(axis.text.x = element_text( hjust = 1, size=a))
CIplot1 <- CIplot1 +theme(axis.text.y = element_text( hjust = 1, size=a))
CIplot1 <-  CIplot1 +  theme(legend.title = element_text(size=16), legend.text = element_text(size=14))

ggsave('plots/WTO/beta1.pdf',CIplot1 , width = 6, height =7.5)

## Figure A17: left panel
nn2 <- c(1, 4, 9, 7, 12, 13); nn3 <- c(1, 4, 9, 7)
upa <- ddd[ c(nn2,(nn2+13), (nn3+13+13)),]
CIplot1 <- ggplot(upa, aes(colour = factor(model)))
CIplot1 <- CIplot1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2)
CIplot1 <- CIplot1 + geom_linerangeh(aes(y = factor(parameter),
                                         xmin = low2, xmax = Upper2), lwd =1, position = position_dodgev(height = 1/2))

CIplot1 <- CIplot1 + geom_pointrangeh(aes(y = factor(parameter)
                                          , x = Mean, xmin = low1, xmax = Upper1), lwd =.45,  position = position_dodgev(height = 1/2),
                                       fill = "WHITE")
CIplot1 <- CIplot1 + theme_bw()
CIplot1 <- CIplot1 + theme(legend.position="bottom")
CIplot1 <- CIplot1 + labs(colour = "Model")
CIplot1 <- CIplot1+ scale_colour_discrete(labels = c("rho-MF", "rho_t-FE",   "rho_t-MF")) + xlab(expression(beta)) + ylab("")+ggtitle("")
CIplot1
a <- 16
CIplot1 <- CIplot1 +theme(axis.text.x = element_text( hjust = 1, size=a))
CIplot1 <- CIplot1 +theme(axis.text.y = element_text( hjust = 1, size=a))
CIplot1 <-  CIplot1 +  theme(legend.title = element_text(size=16), legend.text = element_text(size=14))
ggsave('plots/WTO/beta2.pdf',CIplot1 , width = 6, height = 7.5)

## Figure A17: right panel
upa <- ddd[-c(nn1,(nn1+13),(nn1+13+13), nn2,(nn2+13), (nn3+13+13)),]

CIplot1 <- ggplot(upa, aes(colour = factor(model)))
CIplot1 <- CIplot1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2)

CIplot1 <- CIplot1 + geom_linerangeh(aes(y = factor(parameter),
                                         xmin = low2, xmax = Upper2), lwd =1, position = position_dodgev(height = 1/2))

CIplot1 <- CIplot1 + geom_pointrangeh(aes(y = factor(parameter)
                                          , x = Mean, xmin = low1, xmax = Upper1), lwd =.45,  position = position_dodgev(height = 1/2),
                                      fill = "WHITE")
CIplot1 <- CIplot1 + theme_bw()
CIplot1 <- CIplot1 + theme(legend.position="bottom")
CIplot1 <- CIplot1 + labs(colour = "Model")
CIplot1 <- CIplot1+ scale_colour_discrete(labels =  c("rho-MF", "rho_t-FE",   "rho_t-MF")) + xlab(expression(beta)) + ylab("")+ggtitle("")
a <- 16
CIplot1 <- CIplot1 +theme(axis.text.x = element_text( hjust = 1, size=a))
CIplot1 <- CIplot1 +theme(axis.text.y = element_text( hjust = 1, size=a))
CIplot1 <- CIplot1+  theme(legend.title = element_text(size=16), legend.text = element_text(size=14))
ggsave('plots/WTO/beta3.pdf',CIplot1, width = 6, height = 7.5)


## making graphs of the network
#library(igraph)
#exm <- c(1, 21, 38)
#year <- c(1971, 1992, 2008)
#for (i in 1:3){
#    k <- exm[i]
#    g6 <- graph_from_adjacency_matrix(W3[,,k], weighted=TRUE, mode="upper")
#    E(g6)$width <- E(g6)$weight*10
#    V(g6)$label <- NA
#    V(g6)$size <- 6
#    l <- layout_with_fr(g6)
#    plot(g6,  layout=l)
#    pdf(paste('plots/WTO/Wto',year[i] , ".pdf", sep=""), width=10, height=10)
#    plot(g6, layout=l, edge.arrow.size=0.2, arrow.width=0.2)
#    dev.off()
#}

print(Sys.time()-begin.time) 
sink() # close log  




